Анализ фармакокинетики этанола исследования Yelland-2008

Настройка темы

theme_custom <-
  theme(
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20, margin = margin(t = 15, unit = "pt")),
        axis.title.x = element_text(size = 24),
        axis.title.y = element_text(size = 24),
        panel.background = element_rect(fill = "white", colour = "white"),
        panel.grid = element_line(color = "gray85"),
        legend.text = element_text(size=20),
        legend.title = element_text(size=24),
        plot.title = element_text(size = 30, margin = margin(b = 10, unit = "pt")),
        plot.subtitle = element_text(size = 26, margin = margin(b = 15, unit = "pt"))
  )
standard_theme <- theme_bw() +
  theme(plot.title = element_text(size = 30, hjust = 0.5),
    plot.subtitle = element_text(size = 25, hjust = 0.5),
    plot.caption = element_text(size = 20),
    axis.title = element_text(size = 25),
    axis.text.y = element_text(size = 20),
    axis.text.x = element_text(size = 20),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 22),
    strip.text = element_text(size = 22))
theme_set(standard_theme)
set_flextable_defaults(table.layout = "autofit")

Загрузка и подготовка датасета

data_Yelland <- read.csv("data/raw/Yelland-2008.csv")
data_Yelland <- data_Yelland %>% 
  mutate(across(c(subject_id, occasion), as.character))
str(data_Yelland)

## 'data.frame':    531 obs. of  9 variables:
##  $ study_id     : chr  "Yelland_2008" "Yelland_2008" "Yelland_2008" "Yelland_2008" ...
##  $ subject_id   : chr  "1" "1" "1" "1" ...
##  $ occasion     : chr  "1" "1" "1" "1" ...
##  $ time_min     : int  0 19 30 40 50 59 74 89 109 140 ...
##  $ conc_g_per_L : num  0 0.03 0.4 0.39 0.71 0.66 0.73 0.71 0.67 0.61 ...
##  $ dose_g_per_kg: num  0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 ...
##  $ route        : chr  "PO" "PO" "PO" "PO" ...
##  $ food         : chr  "fed" "fed" "fed" "fed" ...
##  $ site         : chr  "venous" "venous" "venous" "venous" ...

Графики индивидуальных кривых

# current_num <- 1
# for(i in 1:nrow(data_Yelland)) {
#   if(data_Yelland$time_min[i] == 0) {
#     current_num <- 1  # обнуляем счетчик
#   } else {
#     current_num <- current_num + 1
#   }
#   data_Yelland$point_id[i] <- current_num
# }
# 
# rm(i)
# rm(current_num)

# можно сделать то же самое с помощью row_number()
data_Yelland <- data_Yelland %>% 
  mutate(case_name = str_c(subject_id, "_",occasion))

data_Yelland <- data_Yelland %>%
  group_by(case_name) %>%
  mutate(point_id = row_number()) %>% 
  ungroup()

Среднее и стандартное отклонение по всем субъектам (по номеру точки)

Время на точках усреднено

data_Yelland_mean_sd <- data_Yelland %>% 
  group_by(point_id) %>%
  summarise(conc_avg = mean(conc_g_per_L), conc_max = mean(conc_g_per_L) + sd(conc_g_per_L), conc_min = mean(conc_g_per_L) - sd(conc_g_per_L), time_min_avg = mean(time_min))

data_Yelland_mean_sd %>% 
  ggplot() +
  geom_line(aes(x = time_min_avg, y = conc_avg)) +
  geom_errorbar(aes(x = time_min_avg, ymin = conc_min, ymax = conc_max), width = 6, size=0.5) +
  scale_x_continuous(breaks = seq(0, 400, by = 50)) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.25)) +    
  labs(title = "Зависимость концентрации от времени", subtitle = "Mean ± SD", x = "Время, мин", y = "Концентрация, г/л") +
  theme_custom

Расчет основных фармакокинетических параметров

Рассчитаем центральную тенденцию и разброс

pharmacokinetics_data_table <- pharmacokinetics_data %>% 
  mutate(Beta = Beta * 60) %>%  # перевожу в г/л в час
  rename(`β, г/л в час` = Beta,
         `C0, г/л` = C_0,
         `Vd, л/кг` = V_d,
         `AUC0-t, г*мин/л` = AUC,
         `Cmax, г/л` = Cmax,
         `Tmax, мин` = Tmax,
         `Absortion rate, г/л в час` = Abs_rate)
pharmacokinetics_data_table %>%
  select(-subject_id, -occasion) %>% 
  pivot_longer(!case_name,
               names_to = "Показатель",
               values_to = "Значение") %>% 
  group_by(`Показатель`) %>% 
  summarise(across(`Значение`, stat_ph, .names = "{.fn}")) %>% 
  ungroup() %>% 
  flextable() %>%
  theme_box %>% 
  fontsize(size=9) %>% 
  add_header_row(values  = "Описательные статистики фармакокинетических показателей", colwidths = 5)

Описательные статистики фармакокинетических показателей
Показатель	Среднее	Медиана	Стандартное отклонение	Коэффициент вариации (CV), %
AUC0-t, г*мин/л	153.3721	143.1100	34.3078	22.37
Absortion rate, г/л в час	0.0110	0.0103	0.0033	29.74
C0, г/л	0.9852	0.9729	0.1309	13.29
Cmax, г/л	0.7903	0.7900	0.1297	16.41
Tmax, мин	76.1111	74.5000	17.4451	22.92
Vd, л/кг	0.7220	0.7195	0.0898	12.44
β, г/л в час	0.1552	0.1548	0.0186	11.97

Визуализируем выборочные распределения всех параметров из предыдущего пункта

hist_fn <- function(data, col) {
  
  ggplot(data, aes(x = {{col}})) +
  geom_boxplot(
    alpha = 0.5,
  # width = width,
    linewidth = 1.5,
    outlier.size = 2,
  ) +
  geom_histogram(
    fill = "midnightblue",
    alpha = 0.5,
    colour = "black",
    bins = 20
  ) +
    scale_y_continuous(breaks = seq(1, 100, by = 1)) +
    geom_vline(aes(xintercept = mean({{col}}, na.rm = T)),
               color = "darkred",
               linetype = 2,
               linewidth = 2) +
  labs(# title = rlang::englue("{{col}}"),
       y = "Количество, n",
       x = rlang::englue("{{col}}")) +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title = element_text(size = 25),
        axis.text.y = element_text(size = 23),
        axis.text.x = element_text(size = 23),
        legend.title = element_text(size = 23),
        legend.text = element_text(size = 23))
}

pharmacokinetics_data_for_plots <- pharmacokinetics_data_table %>% 
  group_by(subject_id) %>% 
  summarise(across(where(is.numeric), ~mean(.x))) %>% 
  ungroup() %>% 
  arrange(as.numeric(subject_id))

col_list <- pharmacokinetics_data_for_plots %>%
  select(-subject_id) %>%
  names() %>%
  syms()

hist_plots <- purrr::map(col_list,
    ~ hist_fn(pharmacokinetics_data_for_plots, !!.x))

wrap_plots(hist_plots, ncol = 2) +
  plot_annotation(title = "Гистограммы и боксплоты фарм. показателей",
                  theme = standard_theme)

# CVinter и CVintra

ii_in_variability <- function(pharmacokinetics_parametr, subject_id, vis = 3) {
  
  mod <- lm(log({{pharmacokinetics_parametr}}) ~ factor({{subject_id}}))
  
  anova_mod <- anova(mod)
  
  SDb <- sqrt((anova_mod$`Mean Sq`[1] - anova_mod$`Mean Sq`[2])/vis)
  
  CVb <- sqrt(exp(SDb^2)-1) * 100
  
  SDw <- sqrt(anova_mod$`Mean Sq`[2])
  
  CVw <- sqrt(exp(SDw^2)-1) * 100
  
  b_perc <- SDb^2/(SDb^2+SDw^2) * 100
  
  w_perc <- SDw^2/(SDb^2+SDw^2) * 100
  
  return(tibble(
    SDb = round(SDb, 4),
    CVb = round(CVb, 0),
    b_perc = round(b_perc,0),
    SDw = round(SDw, 4),
    CVw = round(CVw, 0),
    w_perc = round(w_perc,0))
    )
}


pharmacokinetics_data_variance <- pharmacokinetics_data_table %>% 
  summarise(across(where(is.numeric), ~ii_in_variability(.x, subject_id))) %>% 
  pivot_longer(everything(),
               names_to = "Показатель",
               values_to = "variability") %>% 
  unnest_wider(variability)

pharmacokinetics_data_table <- pharmacokinetics_data_table %>%
  select(-subject_id, -occasion) %>% 
  pivot_longer(!case_name,
               names_to = "Показатель",
               values_to = "Значение") %>% 
  group_by(`Показатель`) %>% 
  summarise(across(`Значение`, stat_ph, .names = "{.fn}")) %>% 
  ungroup() %>%
  left_join(pharmacokinetics_data_variance)

## Joining with `by = join_by(Показатель)`

pharmacokinetics_data_table %>% 
  flextable() %>%
  theme_box %>% 
  fontsize(size=9) %>% 
  add_header_row(values  = "Описательные статистики фармакокинетических показателей вместе с межсубъектной и внутрисубъектной вариабельностью", colwidths = 11)

Описательные статистики фармакокинетических показателей вместе с межсубъектной и внутрисубъектной вариабельностью
Показатель	Среднее	Медиана	Стандартное отклонение	Коэффициент вариации (CV), %	SDb	CVb	b_perc	SDw	CVw	w_perc
AUC0-t, г*мин/л	153.3721	143.1100	34.3078	22.37	0.1695	17	62	0.1337	13	38
Absortion rate, г/л в час	0.0110	0.0103	0.0033	29.74	0.1721	17	30	0.2603	26	70
C0, г/л	0.9852	0.9729	0.1309	13.29	0.1208	12	85	0.0517	5	15
Cmax, г/л	0.7903	0.7900	0.1297	16.41	0.1323	13	60	0.1070	11	40
Tmax, мин	76.1111	74.5000	17.4451	22.92	0.1421	14	38	0.1811	18	62
Vd, л/кг	0.7220	0.7195	0.0898	12.44	0.1208	12	85	0.0517	5	15
β, г/л в час	0.1552	0.1548	0.0186	11.97	0.1003	10	66	0.0722	7	34

Визуализация отклонений от среднего по субъекту

pharmacokinetics_data <- pharmacokinetics_data %>%
  mutate(Beta = Beta * 60) %>%
  rename(`β, г/л в час` = Beta,
         `C0, г/л` = C_0,
         `Vd, л/кг` = V_d,
         `AUC0-t, г*мин/л` = AUC,
         `Cmax, г/л` = Cmax,
         `Tmax, мин` = Tmax,
         `Absortion rate, г/л в час` = Abs_rate)

β, г/л в час

param_plot_fn <- function(data, col, occasion, subject_id) {
  
  data_upd <- data %>% 
    group_by({{subject_id}}) %>% 
    mutate(subj_mean = mean({{col}}, na.rm = TRUE)) %>% 
    ungroup()
  
  ggplot(data_upd, aes(x = {{col}})) +
  geom_hline(aes(yintercept = {{col}}),
             color = "purple3",
             linewidth = 2) +
  geom_hline(aes(yintercept = subj_mean),
               color = "forestgreen",
               linetype = 2,
               linewidth = 2) +
  labs(title = rlang::englue("{{col}}"),
       y = rlang::englue("{{col}}"),
       x = "") +
   facet_grid(rows = vars(as.numeric({{subject_id}})),
              cols = vars({{occasion}}),
              labeller = labeller(
                .rows = function(x) paste("Номер субъекта:", x),
                .cols = function(x) paste("Номер визита:", x)
      )
) +
    theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 35),
        axis.title = element_text(size = 30),
        axis.text.y = element_text(size = 30),
        axis.text.x = element_text(size = 30),
        strip.text = element_text(size = 24)
        )
}

param_plot_fn(pharmacokinetics_data, `β, г/л в час`, occasion, subject_id)

C0, г/л

param_plot_fn(pharmacokinetics_data, `C0, г/л`, occasion, subject_id)

## Объём распространения, л/кг

param_plot_fn(pharmacokinetics_data, `Vd, л/кг`, occasion, subject_id)

Площадь под кривой 0-t, г*мин/л

param_plot_fn(pharmacokinetics_data, `AUC0-t, г*мин/л`, occasion, subject_id)

Cmax, г/л

param_plot_fn(pharmacokinetics_data, `Cmax, г/л`, occasion, subject_id)

Tmax, мин

param_plot_fn(pharmacokinetics_data, `Tmax, мин`, occasion, subject_id)

Скорость абсорбции, г/л в час

param_plot_fn(pharmacokinetics_data, `Absortion rate, г/л в час`, occasion, subject_id)